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ABSTRACT 



Measurements of primordial non-Gaussianity (/ml) open a new window onto the 
physics of inflation. We describe a fast cubic (bispectrum) estimator of /jvx, using 
a combined analysis of temperature and polarization observations. The speed of our 
estimator allows us to use a sufficient number of Monte Carlo simulations to characterize 
its statistical properties in the presence of real world issues such as instrumental effects, 
partial sky coverage, and foreground contamination. We find that our estimator is 
optimal, where optimality is defined by saturation of the Cramer Rao bound, if noise 
is homogeneous. Our estimator is also computationally efficient, scaling as 0(A 3 / 2 ) 
compared to the 0(N 5 ^ 2 ) scaling of the brute force bispectrum calculation for sky maps 
with N pixels. For Planck this translates into a speed-up by factors of millions, reducing 
the required computing time from thousands of years to just hours and thus making 
Jnl estimation feasible for future surveys. Our estimator in its current form is optimal 
if noise is homogeneous. In future work our fast polarized bispectrum estimator should 
be extended to deal with inhomogeneous noise in an analogous way to how the existing 
fast temperature estimator was generalized. 

Subject headings: cosmic microwave background, early universe, inflation 



1. Introduction 



In the last few decades the advances in the observational cosmology have led the field to its 
"golden age." Cosmologists are beginning to nail down the basic cosmological parameter s, and 
have started ask i ng questions about the nature of the initia l conditions provided by inflation (jGuth 



1981 



Sato 



1981 



Lindd Il982l : lAlbrecht SteinhardtJ Il982l ) , which apart from solving the flatness 



and horizon problem, also gives a mechanism for producing the seed perturbations for structure 
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form ation flGuth Pilll982l ; IStarobinskvlll982l ; lHawking|ll982l ; lBardeen et al.lll983l ; lMukhanov et al 
1992 ). and other testable predictions. 



The main predictions of a canonical inflation model are: (i) spatial flatness of the observable 
universe, (ii) homogeneity and isotropy on large scales of the observable universe, (iii) nearly scale 
invariant and adiabatic primordial density perturbations, and (iv) primordial perturbations to 
be very close to Gaussian. The cosmic microwav e background (CMB ) data from the Wilkinson 
Anisotropy Probe (WMAP) feennett et ahlbood ). both temperature (jfflnshaw et all bood . bood ) 
and polarization (jKogut et al.l 120031 ; iPage et al.l 120061) anisotropies, have provide d hitherto the 



strongest c onstraints on these predictions (jKomatsu et al.l 120031 : iPeiris et al.l 12003 ; ISpergel et al 
20031 . 120061 ) . There is no observational evidence against simple inflation models. 



Non-Gaussia nity from the simplest inflation models that are based on a slowly rolling scalar 
field i s very small dSalopek fc Bondlll99d . Il99ll : iFalk et alJl99ll : bangui et all 19941 : lAcauaviva et al. 
2003 : lMaldacenall2003l ) : however, a very large class of more general models with, e.g., multiple scalar 
fields, features in inflaton potential, non-adiabatic fluctuations, non-c anonical kinetic terms, devi- 
ations from Bunch-Davies vacuum, among others (jBartolo et al.ll2004l . for a review and references 
therein) generates substantially higher amounts of non-Gaussianity. 

The amplitude of non-Gaussianity constrained from the data is often quoted in terms of 
a non-linearity parameter f^L (defined in section 12.11). Many efficent methods for evealuating 
bispectum of CMB tempe rature anisotropies exist (jKomatsu et al.l 120051 ; ICabella et al.l l2006al : 
Smith &; Zaldarriagal 120061 ) . So far, the bispectrum tests of non-Gau ssianity have not det ected 
any sign ificant $nt, in temperature fluctuations mapped by COBE (IKomatsu et al.l 120021) and 
WMAP dKomatsu et ali bood: ISpergel et alibood : ICreminelli et al.l l2006al lbl: ICabella et al.ll2006bl : 
Chen &: Szapudil 120061 ) . Different models of inflation predict different amounts of /nl, starting 
from O(l) to Jnl ~ 100, above which values have been excluded by the WMAP data already. On 
the ot her hand, some authors have claimed non-Gaussian signatures in the WMAP temperature 



2003 



data (IMukherjee fc Wand |200J; lLarson fc Wandeltl |200J; IVielva et all |200J; IChiang et al. 
20041 ). These signatures cannot be characterized by /jvx and are consistent with non-detection of 

Currently the constraints on the /nl come from temperature anisotropy data alone. By also 
having the polarization informat ion in the cosmic microwav e background, one can im prove sen- 



sivity to primordial fluctuations (jBabich Zaldarriaga 



2004 



the experiments have alrady started characterizing pola rization anisotropies (IKovac et al 



Kogut et al. 


2003; 


Page et al. 


2006; 



Yadav fc Wandeltl 1200.^. Alt hough 



2002 



Montroy et al.l l2006). the errors are large in comparison to 



temperature anisotropy. The upcoming experiments such as Planck will characterize polarization 
anisotropy to high accuracy. Are we ready to use future polarization data for testing Gaussianity 
of primordial fluctuations? Do we have a fast estimator which allows us to measure Jnl from the 
combined analysis of temperature and polarization data? 



In this paper we extend the fast cubic estimator of Jnl from the temperature data (jKomatsu et al 
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20051 ) and derive a fast way for measuring primordial non-Gaussianity using the cosmic microwave 



background temperature and polarization maps. We construct a cubic statistics, a cubic combina- 
tion of (appropriately filtered) temperature and polarization maps, which is specifically sensitive 
to the primordial perturbations. This is done by reconstructing a map of primordial perturbations, 
and using that to defi ne our estimator. We also sh ow that the inverse of the covariance matrix for 
the optimal estimator (jBabich &i Zaldarriagal 12004 ) is the same as the product of inverses we get in 
the fast estimator. Our estimator takes only iV 3 / 2 operations in comparison to the full bispectrum 
calculation which takes iV 5 / 2 operations. Here N refers to the total number of pixels. For Planck 
N ~ 5 x 10 7 , and so the full bispectrum analysis is not feasible while ours is. 



2. Primordial Non-Gaussianity 
2.1. A Model 



The harmonic coefficients of the CMB anisotropy a\ r 
to the primordial fluctuation as: 



a £m 



2b 



-± / k l dk r l dr [ $, ro (r) + S em (r) ] h(kr) + n t . 

TT 



T- 1 f d 2 nAT(h)Y* m can be related 

(1) 



where <3?£ m (r) and Si m (r) are the harmonic coefficients of the primordial curvature perturbations 
and the primordial isocurvature perturbations respectively at a given comoving distance r = |r|; 
9xt( r ) is the radiation transfer function of either adiabatic or isocurvature perturbations; X refers 
to either T or E; and ji(kr) is the Bessel function of order I. A beam function bi and the harmonic 
coefficient of noise n^ m are instrumental and observational effects. Eq. ([1]) is written for flat 
background, but can easily be generalized. 

Any non-Gaussianity present in the primordial perturbations <3?£ m (r) or Sg m (r), can get trans- 
fered to the observed CMB i.e. a ; m , via Eq. (fH). Due to the smallness of isocurva ture contribution 



over the curvature perturbation (jPeiris et al 



20031 : iBean et alJl2006l : lTrottall200a ) we will drop the 



isocurv ature contribution f r om Eq. CD) and fu r ther we will use a popular and simple non-Gaussianity 



model (jGangui et al J 1 1994 : 1 Verde et alj|2000l : iKomatsu fc Spergelll200ll ) given by 



Hr) = <S> L (r) + f NL [$l(r)-(<S>l(r))} 



(2) 



where $L( r ) is the linear Gaussian part of the perturbations, and Jml is a non- linear coupling 
constatnt characterizing the amplitude of the non-Gausianity. The bispectrum in this model can 
be written as: 



mhMk 2 Mk 3 ) = 2f NL (2ny6 6 (k 1 + k 2 + k 3 )[P{k 1 )P(k 2 ) + cycl. 



(3) 



The above form of the bispectrum is specific to the model chosen and so in gen eral the con- 
strains on /jvl do not necessarily tell us about non-Gaussianity of other models ()Babich et al. 
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2004k iBartolo et al.l 12004). This model, for instance, fails completely when non-Gaussianity is lo- 



calized in a specific r ange in k space, the case that i s predicted from infla t ion models with f eatures 



in inflaton potential ([Wang fc Kamionkowskil |2000| ; iKomatsu et al.l 12003 ; IChen et al.l 120061 ) . Even 



for the simplest inflation models based on a slowly rolling sc alar field, the b i spectrum of infla - 



ton perturbations yields a non-trivial scale dependence of /tvl (|Falk et al.lll993l ; iMaldacenal 12003 ) 



although the amplitude is too small to detect. On the other hand the bispectrum of curvature 



pertu rbations are exactly Gaussian. This contribution actually follows Eq. 



20051 ). Curvaton models also yield the bispectrum in the form given by Eq. 



(Lvth & Bod 


riguez 


(Lvth et al. 


2003Y 



2.2. Reconstructing Primordial Perturbations using Temperature and Polarization 

Anisotropy 



Reconstruction of primordial perturbations from the cosmological data allows us to be more 
sensitive to primordial non-Gaussianity, which is important because non-Gaussianity in the cosmic 
microwave background data does not necessarily imply the presence of primordial non-Ga ussianity. 
The method of reconstruction from temperature data is described in (jKomatsu et al.ll2005r) and that 
from the combined analysis of temperature and polarization data is described in (jYadav Wandelt 
20051 ). where we reconstruct the perturbations <I>(r) using an operator Oi(r). These operators are 
given by: 



where 




a 



TT 



r<TE 



(jTE (jEE 



-1 



(4) 



C, 



XY 



( a fm a im) 



2Kb 



TT 



k 2 dkP<s,{k)gxig Y i{k)+Ni 



XY 



(5) 



2b x 

P?(r) = (^im(r)a^) = -L- / k 2 dk P*(k) g xe (k) j e (kr), 



TT 



(6) 



P$(k) is the power spectrum of the primordial curvature perturbations, and N XY = (a^a^) is 
the noise covariance matrix. The primordial perturbation then can be estimated as: 



(7) 



X=T,E 



Figure 1 shows the improvement in reconstruction due to additional information from the CMB 
polarization. 
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Fig. 1. — Reconstructed primordial perturbation map using only temperature information (left), 
and, using both temperature and polarization information (right). Perturbations are reconstructed 
at the surface of last scattering. 

2.3. Fast Cubic Statistics 

We can const ruct a quantity S pr irn , that is analogous to the KSW fast estimator of f^L from 
temperature data (IKomatsu et al.ll2005l ) but generalize it to include polarization data as well. This 
quantity has a simple i nterpretation in terms of the tomographic rec onstruction of the primoridal 
potential described in (IKomatsu et alJl2005l ; lYadav fc Wanderm2005l ), It is the radial integral of a 
cubic combination of the scalar potential reconstructions, the B terms with the analogous A term. 



As in (jKomatsu et al.ll2005l ). we form a cubic statistics given by 



1 



J pnm 



f sky 



r 2 dr 



d 2 hB(h, r)B(h, r)A(h, r) 



(8) 



where 



B(n,r) ^Y.Y.^y^mPf^YUn), 

ip Im 



(9) 



ip Im 



(10) 



c? t = J k 2 dkgxe(k) je(kr), and f s k y is a fraction of sky. Index i and p can either be T or E. We 
find (see Appendix A) that (S pr i m ) reduces to 

(Sprim) = E ^fNLBZ^iC- l yf^-%{C~ 1 tB^ m ^ (H) 



sky ijkpqr2<t 1 <t 2 <h ^ l£ ^ 3 



where A^ 2 £ 3 is: 1 when 4 ^2 / 4i 6 when l\ = £2 = £3, and 2 otherwise; B^^ im is the 
theoretical bispectrum for /tvl = 1 ; an d is given by: 
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KZT = J Wa / r 2 dr[/3l(r)Pl(r)al(r) + 0g (r)^(r)a^(r) + ^(^(r^r)] (12) 
where 



(24 + 1) (24 + 1) (24 + 1) 



*■«■«■ =2 ^ o o o J (13> 



Since {S pr i m ) is proportional to Jnl, an unbiased estimator of /jvx can be written as: 

$ _ Sp r i m /*| ^ \ 

JNL — 7 : : : : — — r^r (14) 

The most time consuming part of the calculation is the harmonic transformation necessary for Eqs. 
([9]) and (|10p . The fast estimator defined above takes only TV 3 / 2 operations times the number of 
sampling points for r, which is of the order 100. Hence this is much faster than the full bispectrum 
analysis, which scales as iV 5 / 2 . N here is the number of pixels. In the next section we will show 
that this fast unbiased estimator is also optimal, by proving equivalence with a known optimal but 
slow estimator. 




Fig. 2. — Fisher predictions for minimum detectable /tvl at the 1-er level. Left panel: ideal 
experiment. Right panel: Planck satellite. Solid lines: temperature and polarization information 
combined. Dashed lines: temperature information only. Dot-dashed line: polarization information 
only. 
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2.4. Optimality of the fast estimator 



In their recent paper Babich and Zaldarriaga (jBabich &: Zaldarriagal 12004 ) have found an 



optimal estimator for /tvl which minimizes the expected \ 2 given by 

2_ ( rjijk,obs f Tjijk,prim\ -l\ijk { r>pqr,obs f R pqr,prim\ 

X - 2^ Z> l^Ws _ JNLB iihi3 ) (Cov y ?r (# W3 - f N LB hh£3 ) . (15) 

ijkpqr t\tilz 

This optimal estimator is 

V V R ijfc,o6s rC!rw— 1 V-> p>pqr,prim 

JNL = ^ : : : — , (to) 

J F} t 3 k > pnm (Cnv- 1 \ -npqr,prim^ v > 

l^ijkpqr 2-ut x l % l2, D lii 2 l 3 l^ ov )ijk,pqrO ill2l3 

where the index ijk and pqr runs over {TTT, TTE, TEE, EEE}, unlike in the fast estimator case 
where ijk and pqr run over all the 8 combinations {TTT, TTE, TET, ETT, TEE, ETE, EET, EEE}. 
In appendix B we show that the inverse of the covariance matrix for the BZ estimator is the same 
as the product of inverses we get in the fast estimator, Eq. HH hence our estimator is optima^] 



3. Results 

To test optimality of our fast estimator we use Eqs. (fl4"j) and ([8]) to measure f^L from 
simulated Gaussian skies. The error bars on /jvx are then derived from Monte Carlo simulations of 
our fast estimator. The simulated errors are compared with the Cramer Rao bound (F -1 / f.^ky) 1 ^ 2 



where f. q k v is a fraction of sky, and F is the Fisher matrix given by (jKomatsu Spergel 



Babich fc Zaldarriagal 12004 ): 



2001 



rp _ \" \ " r,pqr,prim , fi-lsip i ^-y-lyg i ^-l^kr T}ijk,prim /-, 7 x 

* - 1^ 1^ A „ a txi*i3 t° M° u 2 \y )iz a txhh > \ u > 

ijkpqr h<l2<lz 

where A^^g is: 1 when i\ / 4 / ^3, 6 when l\ = £2 = £3, and 2 otherwise. We have used ACDM 
model with £l c = 0.26, fib = 0.04, Q\ = 0.7, h = 0.7, and a constant scalar spectral index n s = 1. 
Since the contribution to the integral in Eqs. [9l and [TOl comes mostly from the decoupling epoch (for 
our simulations r^ ec = 13.61Gpc), our integration limits are r m j n = 13425 Mpc and r max = 13865 
Mpc, with the sampling dr = 15 Mpc. Refining the sampling of the r integral by a factor of 2 
does not change our results significantly. The results are summarized in figure [3J We separately 
explore the effects of excluding foreground contaminated sky regions and noise degradation due to 
the instrument. For definiteness we have used the published Planck noise amplitudes which are 
described in Table 2, though we ignore the effect of the scanning strategy and noise correlations. 
We shall come back to this point in § HI 



By an analogous argument one can show that the temperature-only estimator in Komatsu et al. (2005) is also 
optimal, not slightly suboptimal as stated. 
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Fig. 3. — Testing optimality of our /jvx estimator. In all panels lines show the optimal Cramer Rao 
bounds given by • Symbols show the 1-er errors on /jvx derived from Monte Carlo 

simulations. Upper panels: Monte Carlo errors as a function of f s k y - The star shows the WMAP 
Kp2 mask. Triangles: straight sky cuts excluding regions of low galactic latitudes. The left panel 
shows only the effect of the mask, while the right panel includes homogeneous white noise and 
beam smoothing at the level of the Planck satellite. Lower left: Monte Carlo errors as a function 
of 

imax- Lower right: Incomplete sky coverage causes excess variancde of the low £ modes of the 
polarization bispectra. We show results as a function of £ m i n , below which multipoles have been 
removed from the analysis of polarization bispectra in two cases: 1) stars show the simluated errors 
in the noiseless case; 2) the dashed line and triangles show a case with homogeneous noise similar 
to the Planck satellite. The variance excess due to overweighting of the polarization modes is cleary 
visible in the noiseless case. This panel demonstrates that excluding the lowest i polarization modes 
from the analysis avoids this variance excess without significant loss of information. 



We find that the low i modes of the polarization bispectrum are contaminated significantly 
when f s ky < 1 . This is illustrated by the lower right panel of Figure EJ where we sum over all 
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the I modes for TTT contribution to the bispectrum (because the temperature bispectrum is not 
contaminated by the sky cut), while varying t m in for the other terms (EEE,TTE,TEE). Our 
results show that one may simply remove the contamination in the polarization bispectrum due 
to the sky cut by removing i < 10 (when f s k y ~ 0.8) without sacrificing sensitivity to /nl- The 
contamination appears to be less when noise is added, as the dominant constraint still comes from 
the temperature data which are insensitive to the contamination due to the sky cut. 

We can explain this behaviour simply in terms of the coupling between spherical harmonic 
modes induced by the sky cut. The key observation is that the low £ polarization spectrum is very 
small in the theoretical model we chose, but if it could be observed it would add information to the 
Inl estimator. Therefore, for a noiseless experiment, an optimal estimator designed for full sky 
work will assign a large weight to the low I polarization modes. Since the low I spectrum rises very 
steeply towards higher I a sky cut creates power leakage from higher I to lower I. This biases the 
low i spectrum significantly and this bias is amplified by the large coefficient the estimator assigns 
to these modes. In a realistic experiment including noise, the low I polarization modes are difficult 
to measure, since noise dominates. Accordingly, the estimator assigns small weights to the low I 
polarization modes and the sky cut has a much less prominent effect. 



3.1. Computational Speed 



Our fast estimator takes only N 3 / 2 operations times the number of sampling points for the 
integral over r, which is of the order 100. He nce this is much faster than the full bispectrum analysis 
as discussed in ( Babich &: Zaldarriagal 12004 ) . which goes as N 5 ^ 2 . N here is the number of pixels. 
For Planck we expect N ~ 5 x 10 7 , so performing 100 simulations using 50 CPUs takes only 10 
hours using our fast estimator, while we estimate it would take approximately 10 3 years to do the 
brute force bispectrum calculation using the same platform. 



Table 1: Planck noise properties assumed for our analysis 





Central frequency (GHZ) 
30 44 70 100 143 217 353 


Angular Resolution [FWHM arcminutes] 
AT/T intensity a [10~ 6 /xK/K] 
AT/T polarization (Q and U) a [10"VK/K] 


33 24 14 10 7.1 5.0 5.0 
2.0 2.7 4.7 2.5 2.2 4.8 14.7 
2.8 3.9 6.7 4.0 4.2 9.8 29.8 



"Average la sensitivity per pixel (a square whose side is the FWHM extent of the beam), in thermodynamic temper- 
ature units, achievable after 2 full sky surveys (14 months). 
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Conclusions 



Starting with the tomographic reconstruction approach ( Komatsu et al.ll2005l ; lYadav Sz Wandelt 



20051 ) we have found a fast, feasible, and optimal estimator of /nl, & parameter characterizing the 
amplitude of primordial non-Gaussianity, based on three-point correlations in the temperature and 
polarization anisotropies of the cosmic microwave background. Using the example of the Planck 
mission our e stimator is faster by factors of order 10 6 than the estimator described by Babich and 
Zaldarriaga (iBabich &: Zaldarriagal 12004 ) . and yet provides essentially identical error bars. 



The speed of our estimator allows us to study its statistical properties using Monte Carlo 
simulations. We have explored the effects of instrument noise (assuming homogeneous noise), 
finite resolution, as well as sky cut. We conclude that our fast estimator is robust to these effects 
and extracts information optimally when compared to the Cramer Rao bound, in the limit of 
homogeneous noise. 

We have uncovered a potential systematic effect that is important for instruments measuring 
polarization with extremely high signal-to-noise on large scales. The inevitable removal of contami- 
nated portions of the sky causes any estimator based on the pseudo-bispectrum to be contaminated 
by mode-to-mode couplings at low I. We have demonstrated that by simply excluding low I polar- 
ization modes from the analysis removes this systematic error with negligible information loss. 

It has been shown that inhomogeneous noise causes cubic estimators based on the psudo- 
bispectrum with a flat weighting to be signi ficantly suboptimal (IKomatsu et al.ll2003l ). A partial 
solution to this problem has been found by Creminelli et al. ( 2006a IbT). where a linear piece has 
been added to the estimator in addition to the cubic piece. It should be straightforward to apply 
their method to our estimator. 

Finally, our reconstruction approach may be extended t o find fast estimators for higher order 
statistics, for example trispectrum based estimators of /jvl (IKogo Komatsul 120061 ) . This is the 
subject of ongoing work. 
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5. Derivation of Fast Cubic Estimator 

We derive an expectation value of the cubic statistics given by Eq. [8] 

{Sprim) = T~ [ r 2 dr [ d 2 h{B(h,r)B(h,r)A(n,r)) (18) 



f sky 

using the form of B and A as given by Eq. [U] and [TU] 



(W = 7 L E E E (^(^i^W 

^ ijkpqr m\m2ra^ 

j d 2 hY hmi (n)Y l2m2 (n)Y hm3 (n) (19) 

which simplifies to 

(Sprim) = 7^ E E ^(C-XiC-^liC- 1 )^ [ T 2 dr0l i (r)Pl{r)al i {r)Bf l2l3 (20) 

J sky ■ ■-,„« n J 



S ^ ijkpqr titziz 



where 



and 



hihb _ ^m^m^) * 4 .3 _ (21) 



R ijfc ( £l ^ ^ 3 I R^ fc f99^l 

^1^3 ~ Z> I m m m J i3 «if2« 3 m 1 m 2 m 3 
mim,2»ri3 \ / 

is the angular bispectrum, and -B^ 2 £ 3mim2m3 = ( a \ imi a i 2m 2 a hm 3 )' ls tlie CMB bispectrum and can 
be averaged as above due to isotropy. In deriving (S pr i m ) we have also used: 

d 2 nY eimi (h)Y hm2 (n)Yi sm3 (h) = I he2 i 3 ( ^ ^ ^ ) (23) 

\ mi m 2 ^3 / 

The theoretical primordial angular bispectrum can be written as : 

BfXi 3 = fNLtifiSr* (24) 

where 

B pqr,pHm = 2IfM J r 2 dr[Pl{r)(3l{r)al{r) + (r)/f x (r)aj, (r) + /£(r)0S,(r)a&(r)] (25) 
Using the above form of the theoretical bispectrum, (S pr i m ) further simplifies to 

(Sprim) = £" E E B « im (-)( C " 1 )t( C_1 )i 9 (^ 1 )i r ^L 3 - ( 26 ) 
* ijkpqr £i£2^3 

Now since Z^ 2fe = G^^^ A£l ^ ' wnere ^lihtz, which is 1 when l x ^ t 2 ^ £ 3 , 6 when 
£i = £2 = ^3, and 2 otherwise. 

= 77- E E ^/^^^(Ojcc-^cc- 1 )^^ ( 27 ) 

ijkpqr l x <l 2 <h 444 
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Proof of Covariance matrix 



In this appendix we prove the equivalence between the optimal estimator given by (jBabich &; Zaldarriaga 



2004 ) and our fast estimator Eqn. ([26]) 



V 6 

ijkpqri^h a/3 ^2^3 



This is analogous to the temperature-only case (jKomatsu et al.ll2005l ). On the left hand side ijk and 



pgr run over all the 8 possible ordered combinations {TTT, TTE, TET, ETT, TEE, ETE, EET, EEE}, 
while for the estimator on the right hand side a and (3 run only over the four unordered combinations 
{TTT, TTE, TEE, EEE}. 

We prove the equivalence between the optimal estimator and the fast estimator for only one 
combination of l\, £2, and £3, as the proof is same for all the combinations. The covariance matrix 
Cov is obtained in terms of Cj T , Cf E , and Cj E (as in equation 7 in Babich and Zaldarriaga 
(2004)) by applying Wick's theorem, 

CoVj'fc r = (J ir (jil(J k P _|_ (jiq(jjr(jkp _|_ QirQjpfjkq _|_ (jipQjrfjkq _|_ (jiqgjpgkr _|_ (jipQjqQkr ^9) 

The covariance matrix above has the form Cov Q( #, where a, (3, run over {TTT, TTE, TEE, EEE}, 
which after simplification gives: 



/ qTTq,ttqTt 3(7 tt (j TT q t e 3C tt C te C te qTEqTEqTe \ 

3C TT C TT C TE 3C EE C TT C TT +6C TE C TE C TT 3C TE C TE C TE + 6C TT C EE C TE 3C TE C TE C EE 

3C TT C TE C TE 3C TE C TE C TE + 6C TT C EE C TE 6C EE C TE C TE + 3C TT C EE C EE 3C TE C EE C EE 

qTEqTEqTe 3(7 tb (J te C ee 3C te C ee C ee C ee C ee C ee 

1 



where C XY = (C r ) XY are the XY elements of the matrix ( ^ TE n ^E I ■ The above matrix 



x C TE c 1 

is nothing but ^(C~ l y p (C~ l y q (C~ 1 ) kr , after the additional permutations with fixed a and (3 are 
summed up. For example, for a = TTE the ijk index runs over the set {TTE, TET, ETT}. This 
completes the proof. 



